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SUMMARY 

A compressible, unsteady, full Navler-Stokes, finite difference code has 
been developed for modeling transonic flow through two-dimensional, oscillating 
cascades. The procedure Introduces a deforming grid technique to capture the 
motion of the airfoils. Results using a deforming grid are presented for both 
Isolated and cascaded airfoils. The load histories and unsteady pressure dis- 
tributions are predicted for the NACA 64A010 Isolated airfoil and compared with 
existing experimental data. Results show that the deforming grid technique can 
be used to successfully predict the unsteady flow properties around an oscil- 
lating airfoil. The deforming grid technique has been extended for modeling 
unsteady flow In a cascade. The use of a deforming grid simplifies the speci- 
fication of boundary conditions. Unsteady flow solutions similar to the 
Isolated airfoil predictions are found for a NACA 0012 cascade with zero Inter- 
blade phase angle and zero stagger. Experimental data for these cases are not 
available for code validation, but computational results are presented to show 
sample predictions from the code. Applications of the code to typical turbo- 
machinery flow conditions will be presented In future work. 


INTRODUCTION 

The analysis of flow around advanced turboprop airfoil sections requires 
methods capable of modeling unsteady, transonic flow. As the number of blades 
Increase, the cascade effects are expected to become more significant. To 
date, most of the flow codes that model unsteady, transonic cascades are line- 
arized potential solvers. While these codes are fast and more practical for 
load predictions, they are not expected to model the true physics of the flow. 
There Is no experimental data available for propfan sections that can determine 
the extent of the unsteady effects. In an attempt to bridge the gap In under- 
standing and modeling unsteady, transonic flow through cascades, a compressible 
Navler-Stokes code has been developed for such applications. The present code 
Introduces a deforming grid technique to capture the blade motion In the cas- 
cade. The use of a deforming grid Is convenient for treatment of the outer 
boundary conditions since the outer boundary can be fixed In space, while the 
Inner boundary moves with the blade motion. This Is desirable for oscillating 
airfoils In a cascade since the outer boundary position must be known. Sample 
calculations are presented for both Isolated and cascaded airfoils and compared 
to experiment when possible. 
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NOMENCLATURE 


C sonic velocity 

CL lift coefficient 

CM moment coefficient about the leading edge 
CP pressure coefficient 

C«o free-stream speed of sound 

c blade chord length 

e total energy of the fluid per unit volume 
GCL Geometric Conservation Law 

g gap-to-chord ratio 

0 Jacobian of transformation 

k reduced frequency based on semichord, wc/2U ro 

M Mach number 

Re Reynolds number based on chord 

r- upstream-running Rlemann Invariant 

s arclength of a grid line In the n-dlrectlon 

t time normalized by c/C m 

Uoo free-stream velocity 

u.v Cartesian velocities normalized by C® 

V total velocity 

w weighting function for grid deformation 
x.y Cartesian coordinates normalized by chord length 
a angle of attack In degrees 

mean angle of attack 
<r| amplitude of pitching 

y ratio of specific heats 

n normal direction of transformed coordinate system 
E chordwlse direction of transformed coordinate system 
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density 
time variable 

airfoil oscillation frequency 


GOVERNING EQUATIONS 

The code Is an extension of the Isolated airfoil code developed by Sankar 
and Tang (ref. 1). This code solves the two-dimensional, unsteady, Reynolds- 
averaged, compressible Navler-Stokes equations on a body-fitted moving coordi- 
nate system In strong conservation form using an ADI procedure. The cascade 
code uses the same ADI procedure to solve for the Interior of the computational 
domain. The outer boundary conditions are modified to model cascaded airfoils 
and are described later In this paper. A complete description of the formula- 
tion for Isolated airfoils Is Included In reference 1, and only a brief outline 
Is given here. 

The two-dimensional, unsteady Navler-Stokes equations can be written as: 

A A A A A 

q +F C +G =R_+S (1) 

t 5 n S n 


where 


A -1 

q a J {p,pu,pv,e} (2) 

and p Is the fluid density; u and v are the Cartesian components of the 
fluid velocity; e Is the total energy per unit volume. The body-fitted 
(S.ti.t) coordinate system Is related to the Cartesian coordinates using 
the following transformation: 


5 = 5(x,y,t) 
n = n (x,y,t) 

T = t 

The Oacoblan of the transformation Is given by: 

J 3 5 x n y " Vy 3 (x^ - x^) 
and the metrics of the transformation are given by: 


* Jy n 



(3) 


(4) 


(5) 
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Standard central differences were used to compute x^,y 5 ,x n and y n which 
were then used to calculate the metrics. 

The F and G terms In equation (1) are the Invlscld flux terms In the 
5 and n-dl rectlons , respectively. The R and S quantities are the cor- 
responding viscous stress terms'. The viscous terms can be omitted to give 
solutions for the Euler equations. The Beam-Warming, block ADI algorithm Is 
used to solve the governing equations. The viscous terms are treated explicitly 
to reduce the computational time. Artificial dissipation Is added to help the 
stability and follows an Improved method used by Sankar and Tang (ref. 1) to 
reduce excessive smearing of embedded vortices. The discretized form of 
equation ( 1 ) becomes: 


Aq 


n +1 

11 


At 


+ X F n+1 

Vu 


+ 5 e ?; 1 

n 1j 


= 


+ 4 S' 
n 


n n 
G E D 1 j 


(6) 


where 


"n +1 "n +1 "n 

Aq 1j = q 1j ~ q 1j 


(7) 


Is the dissipation term, eg controls the amount of dissipation. At Is 
the time step, and 6 ^ and s n are standard central difference operators. 
The P and G terms are nonlinear and are calculated using a Taylor series 
about the "n" time level. The solution Is second-order accurate In space and 
first-order accurate In time. The Baldwln-Lomax, two-layer algebraic model 
(ref. 2 ) Is used to evaluate the eddy viscosity and the entire flow field Is 
assumed to be turbulent. 


GRID 


The flow solver uses a C-grld generated from Sorenson's (ref. 3) GRAPE 
(GRIds about Airfoils using Poisson's Equations) code. Several test cases 
have been presented by Huff, Wu and Sankar (ref. 4) using the present code and 
C-grlds from the GRAPE code for Isolated airfoils. For cascaded airfoils, the 
grids are generated by specifying an outer boundary shape for the C-grld and 
forcing periodicity on the upper and lower boundaries, as shown In figure 1 . 


Solutions for pitching or plunging Isolated airfoils are found by moving 
the entire grid as a rigid body. This procedure Is valid when the outer bound- 
ary remains In the free stream. For cascade analysis, the outer boundary posi- 
tion must be known to apply the appropriate boundary conditions. A deforming 
grid technique Is used In this analysis to locate the position of the grid as 
a function of time. The inner boundary moves with the airfoil, while the outer 
boundary remains fixed In space. The grid lines connecting the Inner and outer 
boundaries of the C-grld are allowed to deform. This avoids difficult Inter- 
polation of the flow properties along the periodic boundaries. The amount of 
deformation Is a function of the distance away from the surface of the airfoil. 
Define a weighting function, w, to be: 


wU.n) 


s( S , n L- 

s ^ ,n max^ 



(8) 
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where s the arclength of a grid line from the airfoil surface (n = 1) to 
some grid point along \ = constant, and n ma x = the outer boundary grid line. 
From equation (6), we see that w = 1 when s(?,n) = 0 and w = 0 when 
sU,n) = sUtimax)* This function Is used to find the amount of grid deform- 
ation for each new time step: 


n+1 


x 




y 


n+1 

U 


" x 1j + W U (Ax 1j ) r1g 
= y 1j + w 1j (Ay 1j ) r1g 


(9) 


where (Axij) r i g and (Ayij)rlg are the spatial differences If the entire grid 
were rotated as a rigid body. Finally, the amount of deformation of the grid 
over one time step can be defined as: 


(Ax 1J*def = W 1j (Ax 1j^r1g 

( 10 ) 

(Ay 1j J def = w 1j^ Ay 1j^r1g 

The total velocity of the flow at any point on the grid Is equal to the grid 
velocity plus the velocity of the flow with respect to the grid. At the outer 
boundary, the grid Is fixed In space and the grid velocity Is zero. The grid 
velocity changes with the amount of deformation, giving the following Cartesian 
grid velocities: 


(x 


(y 


x _ dx _ (Ax 11 ) de f 
1 j t dr At 

. d y. (Ay 1.1 ) de f 

Ij'x = dx = Ax 


(ID 


With the grid velocities known, the flow solver can predict the unsteady 
flow properties at any desired time. However, as noted by Thomas and Lombard 
(ref. 5) and Shamroth and Glbellng (ref. 6), a moving grid requires careful 
treatment of the metric data to account for the time variation of the Jacobian 
of transformation. The Geometric Conservation Law (GCL) was Introduced by 
Thomas and Lombard (ref. 5) as a way to maintain global conservation of a mov- 
ing grid when transforming from a body-fitted coordinate system to the computa- 
tional plane. A grid that remains fixed In space does not need to account for 
the GCL since the time-dependent terms of the mapping are zero. When the grid 
deforms, the change In the Jacobian of transformation from the previous time 
level must be calculated using the GCL. In the present study, the GCL Is both 
Included and excluded from the calculations to Investigate Its Importance for 
the deforming grids. 


BOUNDARY CONDITIONS 

In order to properly model flow through a cascade, the outer boundary 
conditions must capture the physics of the flow field. The technique used 
here Is similar to the approach used by Chlma (ref. 7). The solution models a 
single blade and assumes periodic boundaries along the upper and lower grid 
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lines. This will be true for cascades with zero Interblade phase angles. 

Based on Shamroth, McDonald and Briley's (ref. 8) experience with modeling 
cascades, the metric data Is also forced to be continuous along the periodic 
boundaries. 

The Inlet conditions are assumed to be at free stream total pressure and 
temperature. The desired flow angle Is also a necessary Input to completely 
specify the Inlet boundary conditions. Slmlllar to Jameson (ref. 9) the 
upstream-running Rlemann Invariant based on total velocity Is calculated and 
extrapolated from the Interior to give a nonreflecting boundary condition: 

R- = [v - 02) 

where V = the total velocity 
C = the sonic velocity 

Although this boundary condition Is only valid for steady flows. It Is 
used here as an approximation for the full unsteady boundary condition by 
assuming that the unsteady terms are negligible. In theory, the nonreflecting 
boundary condition allows the upstream boundary to be closer to the blade with- 
out decreasing the convergence rate. With R - known at the Inlet, the total 
velocity and velocity components are known. The static pressure and density 
are found assuming Isentroplc flow, and the energy Is found from the relation: 

e = P + l p(u 2 ♦ v 2 ) (13) 

This gives all Inlet conditions specified In terms of the conservation vari- 
ables used In equation (2). 

The exit boundary conditions extrapolate p, pu, and pv from the Inte- 
rior and specify free stream static pressure to calculate energy. An exit 
Rlemann Invariant cannot be specified since It would vary across a viscous 
wake. 


The remainder of the boundary conditions are Identical to those used for 
Isolated airfoils In reference 1. This Includes solid wall boundaries on the 
airfoil surface and averaging the flow conditions across the slit aft of the 
airfoil . 


RESULTS AND DISCUSSION 

Results from the sample runs are divided Into two categories: (1) a 

pitching. Isolated airfoil to test the deforming grid technique, (2) steady and 
unsteady predictions for a cascade with zero Interblade phase angles and zero 
stagger. The airfoils used In this Investigation are thicker than typical 
sections found on the propfan. They have been chosen for code validation since 
experimental data for thinner sections does not exist. Preliminary calcula- 
tions have been done for the NACA 16 series airfoils and will be presented In 
a future paper. 
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NACA 64A010 Isolated Airfoil 


The NACA 64A010 airfoil has been chosen for testing the deforming grid 
concept because of the extensive data presented by Davis and Malcom In 
reference 10. This data was taken In the 11 by 11 foot transonic wind tunnel 
at the NASA Ames Research Center. A test case commonly used by other 
researchers Is: M = 0.80, Re = 12.56 million, = -0.21°, a] = ±1.02°, and 

k = 0.202 (reduced frequency based on semichord). This case Is used In the 
present study for comparisons with code predictions. 

Figure 2 shows a fully deformed grid superimposed on a simplified C-grld 
for the NACA 64A010 airfoil. The actual grid used In the predictions Is a 157 
by 58 C-grld, with 97 points wrapped around the airfoil and 60 points In the 
wake. The grid Is Initially generated at the mean angle and then allowed to 
deform so that the maximum deformation occurs at the maximum and minimum pitch- 
ing angles. The plot shows two grids: (1) the Initial, rigid grid generated at 
a = otjj) (solid line)., and (2) the fully deformed grid at the maximum pitching 
angle (dashed line). The grid will also deform In the opposite direction to 
the minimum pitching angle. Notice In figure 2 how the outer boundary remains 
fixed, while the Inner boundary moves with the motion of the airfoil surface. 

The grid connecting the Inner and outer boundaries deforms, as Illustrated by 
the space between the solid and dashed grid lines. 

An average outer boundary distance of ten chord lengths was chosen based 
on the results of a grid parameter study In reference 4 using the NACA 0012 
airfoil. The distance of the first grid line off the surface of the airfoil 
Is 0.00005 chords, which should be adequate for the selected Reynolds number. 

Three separate calculations were studied for this test case. The first 
calculation rotates the grid as a rigid body that follows the airfoil motion. 

The second calculation (deforming, with GCL) uses a deforming grid and Includes 
the Geometric Conservation Law In the calculation of the Jacobian. In order to 
determine the Importance of the GCL for deforming grids, the third calculation 
(deforming, without GCL) uses a deforming grid without conserving the time 
dependent terms In the Jacobian calculation. All three calculations were per- 
formed using five cycles of pitching and were found to be periodic In the fourth 
cycles. The mean angle of attack used In the predictions was adjusted to 
an, = -0.19° Instead of = -0.21° used In the experiment. This was done 
to match the mean lift coefficient. (Since the object of this analysis Is to 
obtain load Information, matching the lift and moment coefficients Is more 
useful than only matching the pressures on the surface with the shock wave. 

The results may still need to be corrected for wind tunnel wall Interference.) 

The lift and moment predictions from all three cases are shown In figures 3 
and 4, respectively. The coefficients from a Fourier transform on the fourth 
cycle are shown In table I, with the results plotted In figures 5 and 6 for 
comparison with experiment. The steady-state pressure distributions for the 
adjusted mean angle are compared with experimental data In figure 7. The real 
and Imaginary components of pressure normalized by the amplitude of oscillation 
are plotted In figures 8 and 9, respectively. 

The results are In reasonably good agreement with experiment and are com- 
parable to the results reported by Chyu and Davis (ref. 13), who also used a 
Navler-Stokes solver. Furthermore, the predictions using a rigid grid are 
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very close to the predictions using a deforming grid, which validates the use 
of a deforming grid for unsteady analysis. There are only minor differences 
between the deforming grid calculations when the 6CL Is Included and excluded. 
The time step used for these calculations was small enough that the error from 
excluding the 6CL Is very small. It Is expected that Including the 6CL In the 
formulation will be more Important for grids with higher deformation and larger 
time steps. Also, the deforming grid technique Is expected to lose accuracy 
as the pitching angle becomes large since the orthogonality and smoothness of 
the grid Is not as good as the numerical grid generated at the mean angle. 

One way to get around this Is to generate another grid after a few degrees of 
pitching and then continue the deformation technique. 


NACA 0012 Cascade 

The deforming grid technique can be extended to predict the flow field for 
an oscillating cascade. No experimental data Is available to validate such cal- 
culations. The solutions presented here will be a first attempt to capture the 
flow behavior for oscillating cascades using a full Navler-Stokes solver. The 
first case to be considered Is a NACA 0012 cascade with g * 3.6, M = 0.80, 

Re = 3.45 million, and = 0.0°. The steady-state solution for pressure 
distribution Is shown In figure 10, along with predictions from other flow 
solvers In references 11 and 12. Caughey and Jameson (ref. 11) show results 
from a transonic potential solver and present solutions using a 80 by 12 grid 
and a 160 by 24 grid. Farrell and Adamczyk (ref. 12) also use a transonic full 
potential solver for the same case using a 98 by 25 grid. The shock strength 
and location predictions from the present analysis are In good agreement with 
the other solutions. The shock Is smeared, but can be made sharper by adding 
more grid lines In the ^-direction. 

Two cases for a gap-to-chord ratio of one are considered next, with 
M = 0.593 and M = 0.660 (Re = 3.21 million, <*„, = 0.0°). These cases repre- 
sent predictions for the cascade with both subsonic and transonic flow. A 157 
by 40 grid was used for these calculations and Is shown In figure 1. The 
unsteady load histories for the airfoils with a-| = ±2.0° about the mean 
angle, M = 0.593, and k = 0.20 are shown In figures 11 and 12. The solutions 
reach a periodic solution by the third cycle. Results for M = 0.660 are not 
plotted since they look very similar. The Fourier Transform coefficients for 
the first harmonics on the fourth cycle of pitching are Included In table I. 
Predictions for the mean, real and Imaginary components of pressure are pre- 
sented In figures 13 to 15, respectively. The mean pressure distributions 
show a weak shock developing when M = 0.660. The Imaginary components of 
pressure are relatively small for both cases. 

All runs were done on the CRAY-XMP computer at NASA Lewis. The steady- 
state solutions required about 2000 Iterations for convergence and took 
3.8xl0 - 5 sec per Iteration per grid point. The unsteady solutions for the 
NACA 64A010 airfoil used 7854 Iterations to complete one cycle of pitching 
using a time step size of 0.0025. Both of the unsteady runs for the NACA 0012 
cascade required a time step size of 0.004 and took 6545 Iterations per cycle 
for M = 0.593 and 5861 Iterations per cycle for M = 0.660. The unsteady 
calculations use 4.3xl0 -5 sec per Iteration per grid point. 
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CONCLUSIONS 


A compressible, unsteady, full Navler-Stokes code has been developed for 
the analysis of oscillating cascades. The present study Introduces a deforming 
grid method to model the airfoil motion and boundary conditions. Results from 
the deforming grid computations are In good agreement with Isolated airfoil 
predictions that use a standard rigid grid. Including the Geometric Conserva- 
tion Law In the Jacobian of transformation has little effect on the solutions 
Investigated In this study. The deforming grid technique has been extended to 
model cascaded airfoils and Is ready for comparison with experiment. To this 
author's knowledge, this Is the first time that a deforming grid has been 
Included In a Navler-Stokes code for oscillating cascade analysis. In addition 
to the valuable aid to the analysis of unsteady flow through cascades, this 
method can be extended to aeroelastlc modeling of nonrlgld airfoils. 
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TABLE I. - FOURIER COEFFICIENTS FOR LIFT AND MOMENT 
TIME HISTORIES (FIRST HARMONICS) 


Case 

a 0 

a l 

b l 

Cl = 3Q + 3] COS((OT 

) +■ b-| Sin (err) 


NACA 64A010 




Rigid grid 

-0.0283 

0.0995 

0.0242 

Deforming grid, with GCL 

-.0283 

.0972 

.0254 

Deforming grid, without GCL 

-.0284 

.0969 

.0255 

Experiment (ref. 6) 

-.0290 

.0862 

.0444 

NACA 0012 




Cascade, M = 0.593 

.0004 

.0708 

-.0323 

Cascade, M = 0.660 

.0006 

.0717 

-.0338 

Cfw| = 3 q + a-| cos (cjt 

) +■ Sin (err) 


NACA 64A010 




Rigid grid 

0.0092 

-0.0352 

0.0015 

Deforming grid, with GCL 

.0092 

-.0344 

.0009 

Deforming grid, without GCL 

.0092 

-.0343 

.0008 

Experiment (ref. 6) 

.0040 

-.0249 

.0004 

NACA 0012 




Cascade, M = 0.593 

.0004 

-.0069 

.0136 

Cascade, M = 0.660 

.0004 

-.0019 

.0146 
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FIGURE 2. - DEFORMING GRID TECHNIQUE FOR AN ISOLATED AIRFOIL. 
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FIGURE 3. - CL VERSUS TIME. NACA 64A010 
AIRFOIL; M = 0.80; d m = -0.19; CL, = 
±1.02; k = 0.202; Re = 12.56X10 6 . 


A 


'RIGID GRID. CM = -0.0257 
DEFORMING GRID (WITH 
. GCL). CM = -0.0250 
DEFORMING GRID (WITHOUT 
„ GCL), CM = -0.0208 



FIGURE 0. - CM VERSUS TIME. NACA 60 AO 10 
AIRFOIL; M = 0.80; a m = -0.19; a i = 
±1.02; k = 0.202; Re = 12.56X10 6 . 




FIGURE 5. - FIRST HARMONICS OF FOURIER 
TRANSFORM FOR LIFT COEFFICIENT ON 
THE FOURTH CYCLE. NACA 60A010 AIR- 
FOIL; M - 0.80; a m = -0.19; a 1 = 
±1.02; k = 0.202; Re = 12.56x10 6 . 


FIGURE 6. - FIRST HARMONICS OF FOURIER 
TRANSFORM FOR MOMENT COEFFICIENT ON 
THE FOURTH CYCLE. NACA 6HA010 AIR- 
FOIL; M = 0.80; a m = -0.19; a, = 
±1.02; k = 0.202; Re = 12.56X10 6 . 
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FIGURE 7 . - MEAN PRESSURE DISTRIBUTION. 
NACA 64A010 AIRFOIL; M = 0.80; 
a m = -0.19; Re = 12.56X10 6 . 


FIGURE 8. - REAL COMPONENT. FIRST HARMONIC 
OF PRESSURE DISTRIBUTION. NACA 6AA010 
AIRFOIL; M = 0.80; a m = -0.19; a 1 = 

±1.02; k = 0.202; Re = 12.56x10 6 . 



FIGURE 9. - IMAGINARY COMPONENT, FIRST HAR- 
MONIC OF PRESSURE DISTRIBUTION. NACA 
69A010 AIRFOIL; M = 0.80; a m = -0.19; 
a 1 = ±1.02; k = 0.202; Re = 12.56X10 6 . 



FIGURE 10. - PRESSURE DISTRIBUTION. NACA 
0012 CASCADE; g = 3.6; M = 0.80; a = 0. 
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FIGURE 11. - CL VERSUS TIME. NACA 0012 CAS- FIGURE 12. - CM VERSUS TIME. NACA 0012 

CADE; M = 0.593; g = 1.0; a m = 0; a 1 * CASCADE; M = 0.593; g = 1.0; a m ■ 0; 

±2.0; k = 0.20; Re = 3.21X10 6 . a, = ±2.0; k = 0.20; Re = 3.21X10 6 . 



CHORD 

FIGURE 13. - MEAN PRESSURE DISTRIBUTION. 
NACA 0012 CASCADE; g = 1.0; a = 0; 

Re = 3.21X10 6 . 



CHORD 


FIGURE 19. - REAL COMPONENT, FIRST HARMONIC 
OF PRESSURE DISTRIBUTION ON UPPER SURFACE. 
NACA 0012 CASCADE; g = 1.0; d m = 0; a 1 = 
±2.0; k = 0.20; Re = 3.21x10 6 . 



FIGURE 15. - IMAGINARY COMPONENT, FIRST 
HARMONIC OF PRESSURE DISTRIBUTION ON 
UPPER SURFACE. NACA 0012 CASCADE; 
g = 1.0; a m = 0; ^ = ±2.0; 
k = 0.20; Re = 3.21X10 6 . 
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